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ABSTRACT 

We investigate the statistics of flux anomalies in gravitationally lensed QSOs as a func- 
tion of dark matter halo properties such as substructure content and halo ellipticity. 
We do this by creating a very large number of simulated lenses with flnite source sizes 
to compare with the data. After analysing these simulations, our conclusions are: 1) 
The flnite size of the source is important. The point source approximation commonly 
used can cause biased results. 2) The widely used Rcusp statistic is sensitive to halo 
ellipticity as well as the lens' substructure content. 3) For compact substructure, we 
flnd new upper bounds on the amount of substructure from the the fact that no simple 
single-galaxy lenses have been observed with a single source having more than four 
well separated images. 4) The frequency of image flux anomalies is largely dependent 
on the total surface mass density in substructures and the size-mass relation for the 
substructures, and not on the range of substructure masses. 5) Substructure models 
with the same size-mass relation produce similar numbers of flux anomalies even when 
their internal mass proflles are different. 6) The lack of high image multiplicity lenses 
puts a limit on a combination of the substructures' size-mass relation, surface density 
and mass. 7) Substructures with shallower mass proflles and/or larger sizes produce 
less extra images. 8) The constraints that we are able to measure here with current 
data are roughly consistent with ACDM Nbody simulations. 

Key words: 






1 INTRODUCTION 

The Cold Dark Matter (CDM) model with a cosmological 
constant (ACDM) has become the standard model of cos- 
mology. This model is in good agreement with a variety of 
observational probes of the large scales distribution of mat- 
ter and galaxies in the Universe and is in general agreement 
with probes of the distribution of mass in galaxy clusters 
and in large galaxies. In the ACDM model, dark matter 
clumps into halos and galaxies form in the halos. On small 
scales, ACDM predicts that dark matter halos exist down to 
very small masses; the exact lower limit depending on the 
properties of the CDM particle and its thermal history. It 
has long been recognized that the number of observed dwarf 
galaxies in the local group of galaxies falls well short of the 
number of predicted halos ( Moore et al.||1999 Klypin et al 



1999 



Diemand et al. 2007b Springel et al. 2008). This is 



referred to as the substructure problem. Either galaxy for- 
mation is highly suppressed in small mass halos or ACDM 
needs to be modified in some way by, for example, chang- 
ing the properties of the dark matter particle or the initial 



conditions for the density fluctuation in the Universe. Warm 
Dark Matter (WDM) is a popular alternative. Whether or 
not these small mass halos exist has been one of the most 
pressing unanswered question in cosmology for a decade. 



Metcalf & Madau (2001) demonstrated that if small- 



scale structure exists in the distribution of dark matter it 
will have a strong effect on the magnifications of quasar im- 
ages in strong gravitational lenses. This effect causes the 
fiux ratio between images to disagree with any lens model 
with a smooth distribution of matter. These cases are call 
anomalous flux ratios. A particular case had been studied by 



Mao & Schneider ( 1998 ) and subsequently it was shown that 



anomalies are common in quasar lenses (Metcalf & Zhao 



2002 Dalai & Kochanek 2002). This work and a number 



of subsequent studies (see 'Z ackrisson k, Riehm||2010 for a 
review of the subject) relied on fitting lens models to indi- 
vidual lens systems. It has not yet been shown clearly what 
can be causing these anomalies and what cannot be causing 
them. 

In a parallel approach, we and others have tried to sim- 



2 Metcalf & Amara 



ulate the lenses directly from cosmological Nbody simula- 
tions to determine if they are consistent with the observed 
frequency of flux anomalies ("B radac et al.|20 05; Amara et al. 
^06 i Macci o et al._2QQ6 , Xu et al |2009t . The first study pre- 
dicted a large number of anomalies, but it may have been 
strongly affected by shot noise. The two more recent and 
higher resolution studies found that the substructure in the 
Nbody simulations is not sufficient to cause the observed flux 
anomalies (also the conclusion of MaoetaL (2004)). This 
is largely because of the small number density of substruc- 
tures near the radii where images form (typically around 
10 kpc in projection). These studies relied on only a few 
projections of a small number of high resolution halos. It 
is possible that these results are a statistical fluke or that 
the observed anomalies are largely caused by dark matter 
objects along the line of sight but not inside the halo of 
the primary lens il Metcalf] (| 2005a|b| ) . Answering the question 
of whether the Nbody simulations have enough small-scale 
structure in them to account for the flux ratio anomalies is 
one of the primary goals of this paper. 

It is very difficult to realistically simulate strong QSO 
lenses from an Nbody simulation. The first, and most impor- 
tant, problem is that shot noise from the discrete particles 
has a strong effect on the image magnifications. Roughly, 
the error in the magnification goes as d/j. ^ fi^ /\/Ws where 
/i is the magnification and Ns is the number of particles 
over which the smoothing is done. Since /i can be large, 100 
or larger in the best cases for detecting substructure, the 
amount of smoothing needed to obtain an accuracy of even 
10% is very large. So much smoothing can even smooth out 
the very substructures one wants to detect. Because of this 



Xu et al. ( 2009 ) replace an Nbody simulation with a simple 



analytic model fit to an Nbody simulation. A second prob- 
lem is that the highest resolution simulations do not contain 
baryons. Baryons have a strong effect on the profile of the 
lens and in some cases dominate the mass within one Ein- 
stein radius. The baryons need to be put in "by hand". A 
third problem is that the extremely high resolution simu- 
lations required provide one, or at best a few, dark matter 
halos. Variations between halos make their lensing proper- 
ties and their tendency to produce anomalies very different. 
It will be demonstrated in this paper that only very limited 
conclustions about the CDM model can be drawn from a 
single simultated lens. 

To avoid these problems, we take a different approach 
in this paper. We produce a large number of analytic lens 
models that are meant to reproduce the population of lenses 
expected in the ACDM model. We then determine the fre- 
quency of flux ratio anomalies in these lenses and compare 
it to the observed frequency. We adjust the properties and 
abundance of the substructures to see what kind of substruc- 
ture is consistent with observations. The allowed statistical 
properties of the substructures are compared with the prop- 
erties of Nbody halos. 

All previous studies, except Amara et al. (2006), have 



also suffered from the problem that the sources are treated 
as infinitely small points. The magnification of individual 
images are calculated by taking derivatives of the gravita- 
tional force at the position of the image. It will be shown in 
this paper, that since the physical size of the quasar radio 
or mid-infrared emission regions are similar to the sizes of 
the substructures of interest the point source magnifications 



are not accurate approximations. We use a new, high speed 
lensing code called GLAMER (Gravitational Lensing with 
Adaptive MEsh Refinement) ('Metcalf"2011) that is the first 
one capable of producing a very large number of simulated 
lenses with finite sources in a reasonable amount of time. 
It does this through an adaptive mesh refinement algorithm 
that will be briefly described in section [ZS] 

In section |2] the models and techniques used to create 
simulated lenses are described. In section [3] the results of 
those simulations are discussed. Ways of comparing the re- 
sults to the available lensing data are presented in section [4] 
The results are compared with the predictions of cosmolog- 
ical Nbody simulations in section |5] A summery and discus- 
sion are given in section l6] 



2 LENS SIMULATIONS 

Our approach in this paper is to produce a large population 
of realistic simulated lenses and then compare their statis- 
tical properties to the observed population of lenses. To do 
this, we must develop a model for the population of gravi- 
tational lens that includes the host, galaxy + dark matter 
halo, and the substructures within the host. We will not con- 
sider the effects of companion galaxies with masses roughly 
equivalent to the primary lens in this paper. 



2.1 Host lens model 

There is significant evidence from lensing and X-ray ob- 
servations that early- type galaxies have a r~^ mass pro- 
files (IGavazzi et al.||2007| [Humphrey & Buote||2010| IChu- 



razov et al.|2010 Fukazawa et al.|2006 ) . In accordance with 
this finding, we model the host lenses as Distorted Singular 
Isothermal Ellipsoids (DSIE). The surface mass density for 
this model is 



K{r,e) 
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where the Einstein radius is 



te = 47r 



^ D,D, 
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and the critical surface density is 



Ds 



AttG DiDi. 



(1) 



(2) 



(3) 



(4) 



(5) 



where Di^ Ds and Dis are the angular size distance to the 
lens, to the source and between the lens and the source re- 
spectively. The first part ([2]) is a Singular Isothermal Ellip- 
soids whose lensing properties have been extensively studied 
(see'Korm ann et al.| ( |1994p for example). The deflection an- 
gle and shear caused by the series in (|3| have been worked 
out by Evans & Witt ( 2003 ) , although with different nota- 
tion. 

The perturbations bn are assumed to be of the same or- 
der as the observed perturbations in the surface brightness 
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profile of of early- type galaxies. Typical values for 63 and 64 
are two or three percent, but accurate statistics are not avail- 
able ("Bende r et al]|1988| [Kormendy et al.||2009| ). We draw 
random values from a Gaussian distribution with variance 
0.005 for 62 and 63 and 0.01 for 64. We take n > 4 terms 
to be zero. In the observations, 64 is usually defined with 
the orientation of this mode fixed to the same axis as the 
axis of the elliptical component to define the "diskyness" or 
"boxyness" of the galaxy. Since the alignment has important 
effects on the lensing properties, we relax this requirement 
somewhat and allow 03, 4 to vary from the position angle 
of the elliptical component. The misalignment is normally 
distributed with variance 3 degrees. 

We also include background shear and convergence in 



the model. Dalai & Watson (2005) calculated the expected 



distribution of 7 and k. in an Nbody simulation at poten- 
tial lenses. They found that n and I7I are both roughly log- 
normally distributed with a variance of ^ 0.03. We assume 
this distribution in our model. Analytic estimates by |Kee-| 
ton et al. (1997) are in agreement with this result, as are 
observations ("Koopmans et a H2006| ) . 

The model described above is what will be called the 
"standard" host model. To test how sensitive magnification 
anomalies are to the host model, we perform a series of tests 
where the distortions to the lens are increased. For the "ex- 
tra distorted model" , we triple the variance in the distortion 
modes and decouple their orientation from the orientation 
of the elliptical component. For the "extra shear model", we 
triple the variance in the background shear and convergence. 



2.1.1 Distributions of host properties 

Calculating the expected distribution of the lenses' red- 
shifts, velocity dispersions and ellipticities requires knowing 
not only the source luminosity and redshift distributions of 
lenses and sources, but also the many selection effects that 
might be important. The sample of lenses we wish to com- 
pare our results with were discovered in many different ways 
and do not have a uniform, well defined selection criterion. 
Instead of trying to model these biases, we use the distribu- 
tions of already known lenses when possible. 

For the lens and sources redshifts, we use the observed 
values for the Castles lenseqil There are 60 lenses with mea- 
sured source and lens redshift pairs. We draw randomly from 
these sets of redshifts. The lenses discussed in section [4] are 
a subsample of these. 

To get a sample of host velocity dispersions, a^ we use 
the velocity dispersions from the SLACS lenses ( [Koopmans 
et al.|200'6 ). This sample of 61 lenses is used to make a cumu- 
lative distribution of a. The discrete distribution is linearly 
interpolated to get a continuous cumulative distribution and 
then this is randomly sampled from. In the SLACS sample, 
the measured velocity dispersion of stars and the velocity 
dispersion of the best-fit SIE models have statistically in- 
distinguishable distributions. We choose to use the best-fit 
SIE velocity dispersions. These values range from 160 to 
396 km s"^ 

The axis ratios, /, are sampled independently from the 
SLACS lenses in the same way as the velocity dispersions. 



No possible correlations between the internal structure of 
the lenses and their redshift are reproduced in this sampling. 
The average of this distribution is / = 0.75, the standard de- 
viation 0.14 and the range is 0.37 < / < 0.98. The SLACS 
lenses are at relatively low redshift because of their selec- 
tion criterion, but observations indicate that the internal 
structure of early-type galaxies do not evolve significantly 
between 2; = 1 and (Thomas et al. 2005). 

We consider only four image quasar lenses in this paper, 
while the SLACS lenses include two image lenses. The asym- 
metry of the lens changes the area enclosed in the tangential 
caustic and thus a sample of four image lenses will tend to 
have more asymmetric lenses than a sample that includes all 
multiple image cases. To correct for this bias, we calculate 
the ratio of the area within the tangential caustic to the area 
within the radial caustic (or "cut" in the case of a DSIE). 
The number of sources used for the lens is then proportional 
to this ratio. More circular galaxies will have less lenses in 
the final sample. This corrects for the bias in the SLACS 
lenses relative to the four image quasars. From to ^ 100, 
source positions are used for each lens model. This method 
of using a variable number of sources per lens is something of 
a compromise; ideally one would have a population of lenses 
that reflected the biases and one source per lens, but to do 
this the caustic structure of each lens would need to be cal- 
culated and then many of those with a small cross-sections 
for producing four images would be discarded. This would be 
computationally inefficient. A small number of sources per 
lens means that the population of high cross-section lenses 
will be better sampled, but if the average number of sources 
per lens is set too low all the lenses with small cross-sections 
will have zero sources. We have set the number of sources 
per lens so that lenses with zero sources are rare (^ 1%). 

For each lens model the source centers are chosen to 
randomly cover a region that encloses the region within the 
tangential caustic. Some of these source positions give rise to 
less than four images (when the source intersects the caus- 
tic or is completely outside the caustic) and some give rise 
to more than four images (when caustics structure is more 
complicated). The cases with less than four images are dis- 
carded in the analysis that follow. 



2.2 Substructure model 

We wish to construct a substructure model that reflects the 
expectations we have from Nbody simulation, but is rela- 
tively simple and has a small number of parameters that 
can be varied to measure the agreement or disagreement 
with ACDM. 

Simulations show that the mass fraction in substruc- 
ture within a projected radius increases roughly linearly 



with projected radius (Springel et al.||2008| [Diemand et al. 



2007a|b ). With a SIE mass model, this implies that the sur- 
face mass density of substructure is constant at least near 
the Einstein radius and interior to it. This will be assumed 
in all cases. 

The mass function of subhalos in Nbody simulations is 
found to be a power-law 



dn 
dm 



(6) 



http://www.cfa.harvard.edu/castles/ 



where n is the number of substructures in a halo. ' Springel| 
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et al. (2008) found that a ^ 1.9 up to about 1/10 of the 



halo mass without any resolved lower mass limit in halos as 
a whole. Transforming mass function into a projected mass 
function in 2 dimensions is not straightforward because of 
mass segregation in the host halo. The projected substruc- 
ture number density will be denoted rj and the projected 
mass function will be drj/dm. 

It found that the substructures of different masses are 
distributed within host halos in remarkably similar ways ex- 
cept that at each radius the mass function has an upper 



mass cutoff (Springel et al. 2008). If it were not for this 



mass cutoff the projected mass function (surface number 
density) would have the same slope as the the total mass 
function. Instead, the projected mass function will become 
steeper than a = 1.9 above some mass scale. We represent 
this effect in our model crudely with an upper mass cutoff 
that is smaller than the one found for the complete mass 
function 



drj 
dm 



m 




mmin < m < rrima^ 
otherwise 



(7) 



with a ^ 1.9. This is a crude model that could be improved 
on in the future. 

The maximum mass in the mass function must be a 
function of host halo size. A mass scale for the host can be 
defined as the mass within a fixed radius (M oc cr^) or the 
mass within a radius where the average density reaches a 
fixed threshold (M oc a^). The latter is the one commonly 
used to define the mass of a halo in cosmology although the 
virial radius is generally larger than the radii over which one 
would expect the SIE model to hold. However, if the concen- 
tration of the halos does not vary greatly within the range 
of host lenses then the same scaling would be expected in 
the inner regions. Making the maximum substructure mass 
a fixed fraction of the host halo mass results in 



,(cr) = Mrr 



(8) 



The same scaling is assumed for the minimum mass. Mmin 
is used as an adjustable parameter to change the mass scale 
and test the data's consistency with a mass cutoff as would 
be expected in many alternatives theories to CDM. The nor- 
malizing halo is fixes to a* = 200 km s~^. 

The normalization of the mass function ([7| needs to be 
set. To agree with Nbody simulations, the fraction of mass 
in substructure at a fixed fraction of the virial radius should 
be the same in all halos. Since i^host oc a and (|8| makes 
average mass scale like a^ the normalization must scale like 
a~^ . Explicitly the result is 



drj 
dm 



= '''(v) 



(i-«) 



(9) 



The parameter 77* is then the total surface number density 
of substructures is a host with a = a* and is not a function 
of projected radius. 

Although the mass fraction in substructure at a fixed 
fraction of the halo radius is the same for all lenses, the 
same is not true at the Einstein radius. Since te oc cr^, the 
total surface density at te is independent of a for lenses 
and sources at the same redshift, which makes the mass 
fraction scale as a^ at this radius. As a result, we might 
expect substructure to be more important for larger lenses. 



The internal structure of the substructures is, for sim- 
plicity, a simple power-law with a cutoff radius 



.w 



(2-/3) 







f Rcutirn^y 



r < Rcut{m,a) 
r > Rcut{m,a). 
(10) 



In the classical analytic treatment, the average mass density 
within the tidal radius is proportional to the average mass 
density of the host within the substructure's orbit ( ,Binney| 
& Tremaine 1987). This implies Rcut{m,a) oc m^^*^ if all 
the substructures are at the same distance from the center 
of the host, which we assume. Since the mass density at a 
fixed fraction of the host halo radius is independent of the 
host size, it is expected that this relation is independent of 
the host size: 



Rcut(m,a) = Rn 



M^ 



1/3 



(11) 



Here i?max is a free parameter describing the size of the 
most massive substructures. In a more realistic model, there 
would be a significant scatter in the i?cut-cr-m relation, but 
for our purposes this relation is sufficient. Using the classical 
tidal radius, the three dimensional distance from the center 
of the lens that this cutoff radius corresponds to is 



i^galactic = 4.3 kpC 

'i?max^''VlO' Mo 



X 



1 kpc 



200 km s"^ 
3/2 /-,n9 A/r. \ 1/2 



M^ 



(12) 



Our fiducial model will have i?max = 0.5 kpc and Mmax = 
10^ M© so i?max ^1.5 kpc is a representative distance which 
is, perhaps, optimistically compact. We will vary i?max from 
0.25 kpc to 4.0 kpc. 

It should be noted that the appropriate /^galactic for 
lensing would be significantly smaller than the average 
i^gaiactic for subhalos in general. Most subhalos are at large 
radii ( ^ lOOkpc) because there is so much volume at large 
radii to make up for the lower weighted number density. 
Projecting along the line-of-sight weights the inner regions 
of the halo more. The difference is an order of magnitude or 
more. This means that the substructure that are important 
for lensing will tend to be denser than the overall population. 

In summery, the substructure model has the free param- 
eters a, /3, Mmax, Mmin, i^max, Tj^ and the normalization host 
velocity dispersion a* which we fix at 200 km s~^. However, 
in the simulations described in the following a and Mmax 
are fixed and the remaining parameters are varied. 



2.3 Ray-shooting 

The sources that we wish to use in our simulation have sizes 
of ^ 10 pc and the substructures can have similar sizes. 
Therefore, it is essential that we be able to calculate the 
magnification of finite size sources. This requirement has 
been widely ignored in the literature because it is difficult to 
map the image of a finite source in a short enough amount 
of time to make it possible to create the large number of 
simulated lenses required for this problem. A new code, 
GLAMER, has been developed for this and other applica- 
tions. This code employs a highly optimized adaptive mesh 
refinement scheme which allows the shapes of the images 
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Figure 1. An example of the refined grid for one particular lens 
and source position. The refinements continue below the resolu- 
tion of this plot. The defiection angle is calculated once at the 
center of each grid cell. There are four images of a 10 pc sources 
in this case. At a higher resolution than is visible here, the lower 
right image breaks into two. 
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Figure 2. The critical curve (outer curve) and the caustic for the 
same lens as in figure [l] The substructure mass range is 10^ — 
1O^M0 with a number density of 77* = 0.5 kpc~^ and a size scale 
of i^max = 0.5 kpc. In this case the a = 214 km s"-*^, 
1.34, ^lens = 0.41 and / = 0.8. 



^source 



Figure 3. These diagrams represent the categorization of four 
image QSO lenses. The large dots represent the images while the 
small dot in each panel is the position of the source. The dashed 
curve is the critical curve (curve along which the magnification 
diverges) and the solid curve is the caustic curve (the curve on the 
source plane that bounds the region in which a source has four 
images). The four panels correspond to the four types of lenses. 
They are, clockwise from the upper left, an Einstein cross, a fold 
caustic, a short-axis cusp caustic and a long-axis cusp caustic. 
Generally, when the source is near one of the cusps in the caustic, 
three of the images will be close together. When the source is near 
the caustic but not near a cusp, two of the images will be close 
together. We define the angular separation between images as the 
smallest angle between the lines passing through those images and 
the center of the lens. The image with the two smallest angular 
separations to other images is the central image of the image 
triplet which includes its neighbors. It is possible that the triplet 
is not well defined, but this very seldom happens in practice. 
The singlet image is the remaining image. The triplet's opening 
angle, A^, is the angle between the dotted lines shown in each 
case. When A^ is small, the lens is "cuspy". The categorization 
of observed lenses into long-axis and short-axis can be made by 
comparing the distance from the center of the lens to the singlet 
image, to the distance from the center of the lens to the central 
image of the triplet. If the former is larger, it is a short axis case, 
and if the later is larger it is a long axis case. In our simulations, 
this proves to be a very good discriminator. 



and their area to be calculated rapidly. (Because of surface 
brightness conservation, the area of a uniform brightness 
image is proportional to its magnification.) This allows us 
to make millions of mock lenses with a finite size source in 
a relatively short amount of time. Figure [l] illustrates how 
the grid is refined to find all the images and their areas. 
Figure [2] shows the critical curve and caustic structure for 
one example lens. For more details on this code, see |Metcalf| 
( [20TH . 

The range of positions in which a substructure will make 
a significant change to the magnification of an image de- 
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pends on the mass of the substructure. To optimize calcula- 
tions, small-mass substructures that are far away from the 
lens are omitted from the calculation while more massive 
substructures further from the lens are included. To accom- 
plish this, a mass dependent cutoff radius from the center of 
the lens is used: 



rmax(cr, m) = 2 rE(cr) + Rcut{m) + 



2mr-B{cr) 

"TTZjcritCmin 



1/3 



(13) 



The first two terms ensure that all substructures within two 
Einstein radii plus the radius of the substructure are in- 
cluded. The third term ensures that any substructure close 
enough to cause a perturbation to the lens that is not well 
approximated as a pure shear will be included. The param- 
eter emin controls how large the variation in the shear across 
the Einstein radius are allowed to be. We set this parame- 



ter to en 



10" 



The contribution from substructures or 



companions outside this range is considered to be part of 
the background shear discussed in section |2.1| as part of the 
host lens model. 

For each lens model (host and substructure), the criti- 
cal curves and caustics are found first. There are sometimes 
multiple, disconnected critical curves. The main tangential 
caustic is found by requiring its critical curve to be the one 
that encompasses the most area while also surrounding the 
center of the lens. The area within the tangential caustic is 
calculated and the number of source positions that will be 
used for that lens is calculated as described in section 12.11 
The sources are required to have their centers inside the tan- 
gential caustic, but they are otherwise randomly distributed. 
Because of the finite source size, some images will be merged 
and this results in less than four images. 

Some lenses have more than the four images that the 
undistorted host model alone would predict. Some of these 
additional images are very small and/or so close to another 
image that they would not be observed as separate images. 
We do a rough initial cut in all cases by merging together 
any images with centroids that are less than 0.1 arcsec apart, 
roughly the resolution of the Hubble Space Telescope (HST). 
Further discussion of additional images is given in the next 
section. 

Table [l] lists the simulation runs that were performed. 
They are in batches of 100,000 lenses with fixed substructure 
parameters. The first five sets of simulations have no sub- 
structure in them and are used to evaluate the importance 
of distortions to the host lens model and establish a base- 
line from which to measure the importance of substructure. 
The the parameters for the remaining twelve simulation were 
chosen to explore the importance of particular substructure 
properties for lensing. Set 2 is taken to be a fiducial model. 
This is a somewhat arbitrary choice, but we do believe that 
it is similar to the predictions of Nbody simulations except 
for the internal profile of the substructures which, as will be 
shown, has relatively little effect on the lensing properties. 
Relative to simulation set 2, set 1 has a higher minimum 
mass (and average mass), set 3 has a lower minimum mass, 
set 4 has a smaller source size, set 5 has more compact sub- 
structure (a smaller i?max), sets 6 and 8 have less compact 
substructure and set 7 has a shallower internal mass profile 
for the substructures. In sets 9 and 10, the upper mass cut- 
off is increased to 10^° M© which is about 10% of the host's 
virial mass. Set 9 has more compact substructures than set 



10. The i?max values are set here so that the size-mass re- 
lation is the same as in sets 7 and 8. For eample, a 10^ M© 
substructure has the same size in sets 9 and 7. The rescaling 
is nessisary because the size-mass relation is normalized at 
the maximum mass in each model which changes between 
these models. In sets 11 and 12, the upper mass cutoff is de- 
creased to 10^ M0. Set 11 has more compact substructures 
than set 12. Again the i?max values are set to preserve the 
mass-size relation between sets 7 and 11, and between sets 
8 and 12. 

The range in surface number density in the simulation 
sets is meant to span the credible range within a CDM-like 
model (Diema nd et al.|2007b||Springel et al.|2008i ). In set 3 
the number density of substructures is much higher for the 
same mass density so because of computer time constraints 
the mass density range for this set does not go as high as 
in the others although the number density goes higher. The 
ranges in 77* are were chosen to cover the realistic range in 
a CDM-like model. 



3 RESULTS 

We create several million simulated lenses and save the im- 
age positions and magnifications. We also store the point 
source magnifications at the centroid of each image and the 
point source magnification for the point in the image that 
is closest to the center of the source. Some of the host lens 
parameters are also stored. In this paper, for ease of com- 
parison, we classify the observed and simulated lenses and 
reduce the position and magnification information to two pa- 
rameters. The parameter A^ is defined in figure l3] A small 
value of A^ indicates the source is near a cusp in the caus- 
tic. Figure[3]also describes what a long- and short-axis lenses 
are. We have found that a good observational way of sorting 
the lenses into these categories is by comparing the angular 
distance between the center of the lens and the singlet im- 
age to the distance between the center of the lens and the 
central image of the triplet. If the former is greater, then the 
lens is a short-axis lens. Otherwise, it is a long-axis lenses. 
The second parameter used to characterize each lens is 

i^_^^±a^ifi±Mi, (14) 

Ml + M2 + Ms 
where "+" is for long-axis lenses and "— " for short-axis 
lenses. The magnifications for the images in the triplet are 
Ml, M2 and /is 5 with M2 being for the central image. The 
original motivation for this parameter was that i?cusp -^ 
asymptotically as a point source approaches a cusp in the 
caustic ( [Schneider fc Weiss|1 992). The i?cusp parameter has 
been widely used because of this model independent predic- 
tion. In practice, i?cusp is not constrained to a very small 
region around zero because of finite source effects and the 
invalidity of the lowest-order expansion of the lensing equa- 
tion around the cusp. And, as will be shown, the distribution 
of i^cusp is not very model independent. 

Figures [4] and |5] shows the distribution of i?cusp and A^ 
for the sample of simulations listed in the captions. It can 
be seen that the simulated lenses occupy a well localized 
regions in these diagrams when no substructure is present. 
Even when substructure is present at the levels investigated, 
the majority of lenses occupy the same regions with a smaller 
number of cases spread out in tails to the distribution. 
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Table 1. Simulation Runs. The top sections show models without substructure, where the source size, distortion and level of external 
shear is varied. The bottom rows show simulation runs with substructure, a = 1.9 in all cases. 



set 



host model 



Mrr 



(Mo) Mrnin (Mo) i^niax (kpc) ^ T]^ {^PC''^) Rs. 



(pc) number of simulations 



standard 

standard 

extra distorted 

extra shear 

no distortion or shear 



10 
1 

10 
10 
10 



100,000 
100,000 
100,000 
100,000 
100,000 



1 

2 
3 

4 
5 
6 

7 
8 
9 

10 
11 
12 



standard 
standard 
standard 
standard 
standard 
standard 
standard 
standard 
standard 
standard 
standard 
standard 



109 
109 
109 
109 
109 

109 
109 
109 

lOio 
lOio 

108 
108 



10« 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 



0.5 
0.5 
0.5 
0.5 
0.25 
1.0 
0.5 
4.0 
1.1 
8.6 
0.23 
1.8 



1 
1 
1 
1 
1 
1 

0.5 
1 
1 
1 
1 
1 



0.013-0.13 
0.013-0.40 
0.013-0.60 
0.013-0.40 
0.013-0.40 
0.013-0.40 
0.013-0.40 
0.013-0.40 
0.013-0.41 
0.013-0.41 
0.013-0.49 
0.013-0.49 



10 
10 
10 
1 

10 
10 
10 
10 
10 
10 
10 
10 



10^ per ?7* 
10^ per ?7* 
10^ per ?7* 
per ?7* 
per ?7* 
per ?7* 
per ?7* 
per ?7* 
per ?7* 
per ?7* 
per ?7* 
per ?7* 



10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 
10^ 



1.2 X 10^ 
3.0 X 10^ 
4.5 X 10^ 
2.0 X 10^ 
3.0 X 10^ 
3.0 X 10^ 
3.0 X 10^ 
3.0 X 10^ 
3.0 X 10^ 
3.0 X 10^ 
3.0 X 10^ 
3.0 X 10^ 



10' r 

10° r 

10"' r 

10"' r 

10"' r 



10' r 

10° r 

10"' r 

10"' F 

10 
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Figure 6. The range in the fractional error, (/ipoint — Mext)/Mext- 
made by using the point source magnification instead of a finite 
size source. 90% of the simulations, in running bins of 5000, fall 
below these curves. The top panel is for a source with a radius of 
10 pc and the bottom panel is for a radius of 1 pc. The dashed 
curves are for no substructure (77* = 0) and the solid curves are 
for ?7* = 0.2 kpc~^ all from simulation set 2. All the images of 
all the four-image systems are used. The black curves are for the 
point magnification calculated at the centroid of the image not in- 
cluding the images that were merged by the 0.1 arcsec merger re- 
quirement. The blue curves are the same but including the merged 
cases. The red curves are for the point source magnification cal- 
culated at the grid point in the image that is closest to the center 
of the source. The errors for the small source are typically at the 
1% level over a wide range of magnifications, when compared to 
nearest point estimates, showing good convergence on the level of 
the numerical noise from the GLAMER ray-tracing code. 



Figure 7. The ratio of the point source magnification to the finite 
size magnification. The shaded regions show where 90% of the 
simulations in running bins of 5000 are, 5% above and 5% below. 
The simulations and color scheme are the same as in figure [6] As 
in figure |6] the upper panel is for a 10 pc source and the lower 
is for a 1 pc source. We see that depending on the method used 
magnification estimates using points will can lead to both random 
errors and biases (due to the asymmetry of the shaded region) as 
compared to the extended source calculation, which is closer to 
the observables. 



Figure |4] shows how important the ellipticity of the host 
lens is to the distribution of i?cusp values. Distortions to the 
SIE model and background shear do broaden the distribu- 
tion, but ellipticity has a particularly strong effect. If only 
low ellipticity lenses are considered the i^cusp values are re- 
stricted to a much narrower band. The sample of lenses is 
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Figure 4. The distributions of A^ and i^cusp that show the importance of distortion to the elhptical lens model and the 
ellipticity distribution of the lenses. The blue regions show histograms. The left column shows the long axis lenses, and 
the left column shows the short axis lenses. The top row is for simulation set "standard" with i^source = 10 pc which has 
random distortions and the full range of ellipticities. The second row is for the simulation set "no distortions or shear" 
which are pure elliptical models with full range of ellipticities. The third row is the same as the first, but with all the models 
with axis-ratio / < 0.7 removed. The fourth row is the same as the second, but with the same axis-ratio cut. The radio and 
infrared observations are shown as red stars. It can be clearly seen that the distribution of i^cusp is highly dependent on 
the distribution of lens ellipticities and that most of the observed i^cusp values are not exceptionally high if the full range 
of ellipticities is considered. The horizontal lines show where 95% of the cases are above and 95% of the cases are below in 
bins of 3000 simulations. The observed lenses are shown in red and discussed in section 0] 



biased toward high ellipticities relative to the general popu- 
lation of lenses because the cross-section for producing four 
images (the area within the tangential caustic) is increases 
with increasing ellipticity. At the same time, Nbody simu- 
lations might be biased toward low ellipticity since gener- 
ally only dynamically well relaxed systems are chosen for 
very high resolution simulations. This can explain some of 
the discrepancies between simulations and observations that 
have been reported ( [Metcalf &: Amara|2010| |Xu et al.|2009 



Maccio et al.||2006 ). This will be further discussed in sec- 
tion|5] 

Figure |5] is similar to figure |4] but the effect of sub- 
structure on the A^-i?cusp distribution is illustrated. An ad- 
ditional 10% error on each image's flux is added to conser- 
vatively account for typical obervational uncertainties. Sub- 
structure has the effect of producing a papulation of extreme 
outliers in this distribution. 

Figure [6] shows the fractional error made in the mag- 
nifications when the point source magnification is used. It 
can be seen there that the fractional error is small for mag- 



nifications less than around 5. This is confirmation that the 
numerical errors made by the ray-tracing code are small. 
At higher magnifications, larger errors are made when the 
source is 10 pc. This is not a numerical effect. It can also be 
seen in figure [6] that substructure causes the errors made by 
using the point magnification to increase when the source 
size is 10 pc, but less so when the source size is 1 pc. This 
is in agreement with expectations because the source size of 
10 pc is closer to the characteristic scale of the substructures. 
Figure [T] shows the ratio between the point source mag- 
nifications and the finite source magnifications. Again, it can 
be seen that numerical errors are not playing a large part. 
It is evident that the point source magnifications are not 
evenly distributed around the finite source magnifications. 
Centroid point source magnifications tend to overestimate 
the real magnification; in some cases by a large factor. This 
is the magnification that would be calculated when fitting a 
lens model to an observed lens. In the simulation, the cen- 
troid is calculated by doing a flux weighted average over 
the pixels on the simulation grid. The nearest point mag- 
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Figure 5. The distribution of A^ and i^cusp for four of the simulations in set 2. The blue regions show histograms. The left 
column shows the long axis lenses, and the right column shows the short axis lenses. The number density of substructures 
in each row from top to bottom are 77* = 0, 0.09, 0.16 and 0.27 kpc""^. Random noise of 10% has been added to represent 
observational errors. The horizontal lines show where 95% of the cases are above and 95% of the cases are below in bins of 
3000 simulations. These are not exactly the same bins as that are used in calculate the outliers discussed in sect ion [3. 1| but 
are similar. The radio and infrarad data are shown as red stars. 



nification is much less biased and in the opposite direction; 
the magnification is underestimated. In other lensing simu- 
lations, the source position is often fixed and the images are 
found by an iterative minimization algorithm. This would 
give essentially the same result as our nearest point magni- 
fication. Both effects are much smaller for a smaller source 
size, as they should be. 



Many images were merged because their centroids were 
within 0.1 arcsec. In these cases, it makes no sense to take 
the closest point magnification since the closest point is 
not unique. Unsurprisingly, the magnification at the cen- 
troid point is an even worse approximation in these cases, 
as can be seen in figures |6] and |7| In exceptional cases, the 
centroid might not even be in one of the images that are 
merged. As expected, these cases only arise when substruc- 
ture is present. 



Figures |6] and |7| should give one pause before using the 
point source approximation for the magnification in any sub- 
structure lensing study or when interpreting the results of 
any studies that use this approximation. 



3.1 frequency of A^ 



RcMsp outliers 



To determine how often it would be expected for a lens 
to have A^ and i^cusp values that are inconsistent with a 
smooth lens model we define a region around the distribution 
in the case where no substructure is present and find how 
many simulated lenses lie outside this region when substruc- 
ture is added. We define this region by taking bins in A^ that 
contain 2000 simulations taking the long-axis and short-axis 
cases separately. Upper and lower boundaries within each 
bin are set such that 2.5% of the simulations in the bin are 
greater than the upper bound and an equal number are less 
than the lower bound. The bins completely cover the full 
possible range of A^. Without substructure, 5% of a lens lie 
outside of this region. The fraction of simulated lenses out- 
side this region when substructure is added will be called 
the fraction of outliers. 

Figure [S] shows the fraction of outliers as a function of 
the substructure surface number density, 77*, for different 
substructure minimum masses (simulation sets 1, 2 and 3). 
A significant fraction of the lenses are found to be outliers. 
The top panels of figures [O] through pT] show the same outlier 
fraction, but as a function of surface mass density. 

It is surprising that in figure [9] the outlier fraction ap- 



10 Metcalf & Amara 



pears dependent only on the total surface mass density and 
not on the lower mass cutoff. One might think that all the 
lensing is being done by the most massive substructures 
and this is why the lower mass cutoff is not important in 
these cases. This does not seem to be the case; from set 1 
(Mmin = lO^Mo) to set 3 (Mmin = lO^M©) the mass den- 
sity in the highest decade of mass (10^ to 10^ M©) drops 
by 60% for the same total surface mass density and yet the 
number of outliers is unchanged. 

Figure ^] shows the importance of compactness and 
internal structure on the number of outliers. The substruc- 
ture mass function is the same for all the models in this 
figure. The slope of the internal density profile, /3, seems 
to have very little effect on the outlier fraction. On the 
other hand, the size of the substructures, or their compact- 
ness, does have a strong influence of the outlier fraction. 
Between i?max = 0.5 kpc and i?max = 4.0 kpc the frac- 
tion decreases significantly. Since the size-mass relation of 
the substructures is related to their galactocentric distance 
through tidal stripping, this sensitivity would provide infor- 
mation on where the substructures are within the lens halo 
or outside of it. 

In figure [TT] the upper substructure mass limit is 
changed to investigate further the insensitivity to mass 
range. It is seen again that for the same mass-size relation 
the fraction of outliers is dependent on the total surface 
mass density and relatively insensitive to the upper mass 
cutoff. The sensitivity to substructure compactness is again 
clearly present. Set 9 with Mmax = 10^*^ M© appears to pro- 
duce slightly less outliers than set 2 with Mmax = 10^ M©. 
This could be because large substructures will sometimes 
displace the image positions and magnifications significantly 
while preserving a low i^cusp value; the cusp in the caustic 
is moved, but its shape remains relatively intact. 

From the upper panels of figures [5] through [TT] it can 
be seen that if the size-mass relation is held fixed the outlier 
fraction is largely a function of the total surface mass density 
in substructures and not the range of substructure masses. 
This conclusion may depend on the function used here {a = 
1.9). Further simulations will be needed to investigate this. 
Changing the size-mass relation so that the substructures 
are less dense does reduce the fraction of anomalies (sets 6 
and 8). 



4 COMPARISON WITH DATA 

To avoid contamination from microlensing by stars in the 
lens galaxy, differential extinction and variability of the 
source on time-scales smaller than the image delay times, 
we compare our simulations only to quad lenses measured 
in the radio and the mid-infrared. Since we have not included 
companion galaxies to the primary lens in our simulations, 
we also remove lenses with nearby galaxies that appear to 
have similar masses to the primary. This removes 1608+656 
and 1004+4112 from the list. There is a very faint dwarf 
galaxy within the Einstein radius of 2045+265 ("McKean" 
et al.|[2007 ), but we will consider this to be a substructure 
and not a companion galaxy because it is small. Lens models 
show that this substructure would need to be unnaturally 
elongated to cause the flux anomaly in this system, so there 
is probably another substructure present. The lenses must 
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Figure 8. The fraction of lenses that lie outside the region in 
A^-i?cusp space contains 95% of the lenses when there is no 
substructure (see text for details) as a function of substruc- 
ture number density. In all cases Mmax = 10^ M©. The curves 
are for Mmin = lO^M© (red), Mmin = lO'^M© (green) and 
Mmin = 10^ M0 (blue,). These correspond to sets 1, 2 and 3, re- 
spectively from Table [l] There are 100,000 simulated lenses used 
in calculating each point. 
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Figure 9. The top panel is the same as in figure Is] except it is 
now as a function of the surface mass density in substructures. 
The bottom panel shows the fraction of lenses with more than 4 
images with separations of more than 0.1 arcsec and flux ratios 
of within a factor of 100 (excluding the cases with less than four 
images). The colors are the same as in figurels] The dotted lines in 
the bottom panel show where there is only a 10% and 5% chance 
of a sample of 32 lenses having no cases of more than 4 images as 
in the Castles lens sample. 
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Figure 10. Same as figure [9] but with substructure models with 
different internal structures. Green is set 1 (see table [TJ as in 
figure I9] Red is for denser substructures (set 5), while blue and 
purple are for less dense substructures (sets 6 and 8 respectively). 
Cyan is for substructures with less steep mass profiles, but the 
same sizes as green (set 7). Changing the internal mass profile has 
no discernible effect on the frequency of fiux anomalies, but has 
a significant effect on the frequency of high multiplicity lenses. 



Table 2. Observed lenses used in this analysis 



band 



AO 



Re 



reference 



2045+265 


radio 
radio 
radio 

mid-IR 
radio 
radio 

mid-IR 


35.3° 

79.8° 

108° 

74.9° 

101.5° 

146.3° 

127.5° 


0.501 
0.254 
0.417 
0.203 
0.220 
0.357 
-0.043 


Fassnacht et al. (1 


999f 


0712+472 


Jackson et al. ( 19£ 


1555+375 


Marlow et al. ( 1999 ) 


1422+231 
0414+053 
2237+030 


Chibaet al. (2005) 
Katz et al. (1997) 
Falco et al.,(,1996^ 


1115+080 


Chiba et al. 


2005 


' 



also have a detected lens galaxy which eliminates 0134-0931 
and 0128+437. Table [2] lists the lenses used and their i^cusp 
and A^ values are plotted in figures |4] and [5] 

The most striking thing in figures [4] and |5] is that one 
of the lenses, 2045+265, has significantly higher i?cusp than 
is expected in the absence of substructure, but that all the 
other lenses have A^-i?cusp values that are not particularly 
anomalous. In the bottom two rows of figure |4] it can be seen 
that if only the low ellipticity lenses (axis ratio > 0.7) were 
considered three or four of the observed lenses would have 
anomalous A^-i?cusp values. Since the authors that have 
compared Nbody simulations to the data using i^cusp values 
in the past have used very few simulated lenses and all with 
axis ratios ^0.7 (Amara et al.|2006 Maccio et al.|2006 Xu 



et al.|2009 ) it is now not surprising that they concluded that 
the simulations did not produce enough anomalies 



Figure 11. Same as figure [9] but with substructure models 
meant to explore the importance of the upper mass limit. Green 
is set 1 (see table fl]) as in figure [9] and flO] The purple and 
brown (sets 9 and 10, respectively) have a higher upper mass 
cutoff of Mmax = IO^'^'-'M© and different mass-size relations. The 
red and blue curves (sets 11 and 12) are for a mass cutoff of 
Mmax = 10^ M0. The i^max values are set so that a substructure 
of the same mass will have the same size in sets 1, 9 and 11. And 
the same is true for sets 10 and 12. The compactness clearly has a 
strong effect. The fraction of outliers for the green, red and purple 
are similar indicating that the upper mass cutoff does not have a 
strong effect on it when the substructures are compact. In con- 
trast, the number of high multiplicity lenses clearly is dependent 
on the mass range. 



It should be emphasized that just because the lenses' 
A^-i?cusp values are not anomalous does not mean that 
they do not have anomalous flux ratios. Some of these cases 
clearly cannot be fit by reasonable models without substruc- 
ture when all the image positions and fluxes are taken into 



account (Metcalf & Zhao 20021 Evans & Witt 2003 Shin 



& Evans 2008). With so few observed lenses and only one 
clear anomaly in A^-i?cusp space, it is impossible to make 
any strong conclusion about the aloud properties for sub- 
structure using only the A^-i?cusp distribution. About one 
anomalies out of the seven lenses is about what one would 
expect from studying the top panels of figures |9] and [lO] for 
a substructure surface density of ^ 10^ M© kpc~^. Other 
flux-based constraints are possible and will be investigated 
in future papers. 

We introduce another constraint in the bottom panels 
of figure [9] through [TT] based on the fraction of simulations 
with more than four images. (This does not include the cen- 
tral demagnified image that forms near the center of the lens 
for nonsingular lens mass profiles. In our case, the mass den- 
sity in the center of the lens diverges like S oc r~^, and this 
image never appears; it is infinitely demagnified.) Even af- 



12 Metcalf & Amara 



ter merging images with centroids less than 0.1 arcsec apart, 
there are cases where the substructures cause further spht- 
ting of the images. Of the 32 QSO lenses in the Castles 
(Kochanek et al. 2000) list of lenses with more than four 
images and simple lenses, none have more than 4 images 
of a single source separated by more than 0.1 arcseqj This 
puts a strong constraint on the allowed fraction of lenses that 
have more than four images, />4. The probability of getting 
zero cases of > 4 images in 33, given that the probability of 
getting such a case is p ^ />4, is a binomial distribution. 
There would be less than a 5% chance of this happening in 
the observed sample if />4 is greater than 0.089 and less 
than 10% chance if />4 > 0.069. These are the dotted lines 
in the bottom panels of figures [9] through pT] For equal sur- 
face mass density, more massive substructures cause more 
high image multiplicity lenses. 

The multiplicity constraint does change significantly if 
the resolution cutoff of 0.1 arcsec is changed. There are a 
large number of lenses where the images are merged in some 
cases (up to ^ 20%). With improved resolution or a more 
careful anaylisis of the data, we believe this constaint could 
be made significantly stronger. 

Within the ranges of 77* studied here, the only models 
that are limited by this image multiplicity constraint are 
set 1 (high lower mass cutoff and compact), set 5 (super- 
compact) and set 9 (high upper mass cutoff and compact). 
The constraints are S* < 2.0 x lO'^ M© kpc~^, S* < 1.2 x 
10'^ M© kpc~^ and S* < 1.2 x lO'^ M© kpc~^ respectively. 
The more compact and massive the substructures are the 
more high multiplicity cases are created. This constraint is 
in contrast to the i?cusp constraint which depends only on 
the mass density and compactness. With more lenses this 
constraint could become significantly stronger in the future. 



5 EXPECTATIONS FOR SMALL-SCALE 

STRUCTURE WITHIN THE CDM MODEL 

A good point of comparison between lens simulations and 
Nbody simulations is the fraction of mass in substructure 
within a projected radius of 10 kpc. This is easily measured 
in the simulations and since the Einstein radius is typically 
around 10 kpc, it is close to what is actually constrained by 
the lensing data. In our model, this quantity is given by 

Msuh{R < 10 kpc) _ G{m)r], 



/>10kpc 
J sub 



Mhost(i?< 10 kpc) 
1.08 X 10"^ kpc^ S*, 



-10 kpc, (15) 
(16) 



where the fiducial value a* = 200 km s~ has been used. 
Note that this fraction scales with host mass in our model 
and in the simulation s. 

plOkpc 



Diemand et al. 



and 



( |2007a) give ,t. 



Xu et al, 



Via Lactea simulation, 

0.0025 with a large scatter in the Aquarius simulations. 



0.003 for the 



(2009) give/3^^°^P^ 



■'sub 



^ 0134-0931 might have five optical images, but two of them are 
well within 0.1 arcsec of each other. 1933+503 has 10 radio im- 
ages, but models show that the best explanation is that there are 
3 sources with none of them imaged more than four times ( |Nair| 
|1998| ). 1359+154 does appear to be an honest-to-goodness case 
of a single QSO with 6 images, but the lens is a group of three 
galaxies and thus does not pass our no companions cut. 



These simulations should be resolving substructure to be- 
low 10'^ M0. These translate to S* = 2.8 x 10^ M© kpc~^ 
and 2.3 X 10^ M© kpc~^ respectively. Accounting for the ex- 
tra mass below there resolution and judging from figures |9] 
through \n\ we would expect about a 10% chance of a clear 
outlier in the A^-i?cusp distribution for the high compact- 
ness cases which is consistant with the one out of seven 
observed. For the larger size-mass relation (sets 8, 10 and 
12) the expected fraction is increased by only a few percent 
from the no substructure case, but with only one observed 
outlier, we do not consider this a significant contradition. 

[Amara et al.| ( |2006t , |Macci6 et"aL] ( |2006| ) and |Xu et aT 
( 2009 ) come to the conclusion that the substructure present 
in the simulations is not enough to cause the observed fre- 
quency of i^cusp anomalies. In light of the findings in this 
paper we believe that these conclusions were flawed because 
the full range of host lens ellipticiites was not represented in 
the simulations. Maccio & Miranda (2006) may have used 
too low a substructure mass range (10 — 10^ M©) to cause 
enough anomalies. 

There are a number other complicating factors that 
make comparing observations to the true predictions of 
CDM difficult. For example, the baryons are not accounted 
for in the Nbody simulations. This impacts the predictions 
in several ways. First, the host galaxy needs to be inserted 
by hand into these Nbody simulations for them to be re- 
alistic lenses. The mass fraction decreases with the inclu- 
sion of baryons. Second, the baryons are expected to have 
some effect on the internal structure of the substructures, ei- 
ther expanding or contracting them, which will affect their 
tidal stripping and disruption in the host halo. The resi- 
dent galaxy might also have a significant effect on the sur- 
vival of substructures. As discussed in section [2T] the typ- 
ical galactocentric distance for substructures that are im- 
portant for lensing is significantly smaller than the typical 
distance of substructures in general. The substructure pop- 
ulation probed by lensing is likely to be more compact and 
have a steeper mass function, at least above ^ lO^M©, than 
the general population. This steepening of the mass function 
at high masses has been only crudely accounted for in our 
model by the Mmax cutoff parameter. 

Because we appear to be consistent with the simula- 
tions on the frequency of i^cusp anomalies does not mean 
that some other test, such as fitting each simulated lens to 
a smooth lens model, would not show some inconsistency. 
Modeling the lens puts constraints on the ellipticity. Our 
argument is that i^cusp is not a good test for the existence 
of substructure without further constraints. 



6 CONCLUSIONS &^ DISCUSSION 

We have preformed the largest number of lens simulations 
ever done with finite size sources. This was made possible by 
the new adaptive ray-tracing code GLAMER. We find that 
accounting for the finite size of the source is necessary for 
drawing accurate conclusions from the lensed QSO data. 

We find rough consistency between the ACDM predic- 
tions and observations, i^cusp is found to be a poor discrimi- 
nator between lenses with substructure and without because 
of its sensitivity to the ellipticity of the lens. The distribution 
of ellipticities used in our lens models is based on the ellip- 
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ticities of observed lenses so we do no think the elUpticities 
required to explain the observed i^cusp distribution (accept- 
ing lens 2045+265) are atypical. Other methods for com- 
paring observations to models are likely to be more fruitful. 
And as the data improves more precise comparisons will be 
possible. In addition to the substructure within the primary 
lens, there should be some contribution from intergalactic 
small-scale structure (Metcalf 2005a b) so one should expect 
the limits derived from the data to be somewhat higher than 
the limits derived from Nbody simulations of individual dark 
matter halos. The baryons also clearly play a role in shap- 
ing the lensing properties and they are not fully taken into 
account in the simulations at the necessary resolution. 

We have limited our study here to a substructure mass 
function of the form dN/dm oc m~^ with a = 1.9. This 
seems well motivated by the simulations on small mass- 
scales, but could be steeper on larger mass-scales because 
of tidal stripping and disruption in the central regions of 
the lens. With the a = 1.9 mass function, the smaller mass 
substructures plays a smaller part in causing flux anoma- 
lies because most of the mass resides in larger mass objects. 
This will make it difficult to measure any possible lower mass 
cutoff using monochromatic QSO lensing alone. Fortunately 
there are some other prospects for probing the mass func- 



tion in the future such as spectroscopic gravitation (Mous- 
takas fc Metcalf[ 2003) and Einstein rings (Vegetti & Koop- 
mans|20Q9 ). If the slope of the mass function is steeper than 
a = 1.9, the smaller structures will play a larger role in the 
lensing. 

It is clear that what is really required to make a more 
conclusive measurement of the amount of substructure in 
dark matter halos is more data. With 7 lenses, only limited 
conclusions can be made from a statistical point of view. 
We are also vulnerable to systematic errors. For the kind of 
study done here, more strong lenses measured in the radio 
and/or mid-infrared are needed. Planned large scale imaging 
surveyqj expect to increase the number of lensed QSOs in 
the visible by an order of magnitude so we look forward to 
great improvements in this field. 
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